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ABSTRACT 

In this paper we describe a new algorithm for the long-term numerical 
integration of the two-body problem, in which two particles interact under a 
Newtonian gravitational potential. Although analytical solutions exist in the 
unperturbed and weakly perturbed cases, numerical integration is necessary in 
situations where the perturbation is relatively strong. Kustaanheimo-Stiefel 
(KS) regularization is widely used to remove the singularity in the equations 
of motion, making it possible to integrate orbits having very high eccentricity. 
However, even with KS regularization, long-term integration is difficult, 
simply because the required accuracy is usually very high. We present a new 
time-integration algorithm which has no secular error in either the binding 
energy or the eccentricity, while allowing variable stepsize. The basic approach 
is to take a time-symmetric algorithm, then apply an implicit criterion for the 
stepsize to ensure strict time reversibility. We describe the algorithm in detail 
and present the results of numerical tests involving long-term integration of 
binaries and hierarchical triples. In all cases studied, we found no systematic 
error in either the energy or the angular momentum. We also found that its 
calculation cost does not become higher than those of existing algorithms. By 
contrast, the stabilization technique, which has been widely used in the field of 
coUisional stellar dynamics, conserves energy very well but does not conserve 
angular momentum. 

Subject headings: Numerical simulation, Regularization, Time-symmetric 
integration, Two-body problem 
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1. Introduction 

The long-term numerical integration of the two-body problem is important in several 
areas of astrophysics, particularly in the context of the evolution of star clusters. For 
example, binaries play a crucial role in the secular evolution of globular clusters (see, e.g. 
Hut et al. PASP 1992 and references therein). Binary evolution is also important from 
the observational point of view, since many interesting objects in globular clusters, such as 
X-ray sources, millisecond pulsars, high-velocity stars, and blue stragglers, are believed to 
be the result of binary interactions. In order to study the evolution of binaries in globular 
clusters, self-consistent A'"-body simulation is a most useful tool. In such simulations, we 
integrate numerically the orbits of both the stars and the binaries in the cluster. 

There are, however, technical difficulties in obtaining useful results from numerical 
experiments. If we want to study the evolution of globular clusters, we have to follow the 
orbits of tightly-bound binaries under weak perturbations for, say, 10^° years, corresponding 
to several cluster half-mass relaxation times. Since the period of a binary is typically 
less than one year, this time interval corresponds to more than 10^° binary orbits, or 
~ 10^^ — 10^"^ numerical integration steps. The computational requirements of such an 
undertaking would be prohibitive; in addition, both truncation and round-off error would 
be unacceptably large. Of course, there are several techniques used in order to reduce the 
number of steps. First of all, in practice, we deal with this difficulty by regarding binaries as 
"unperturbed" if the relative perturbation is less than, say, 10~^ or so. By this treatment, 
the relative calculation cost of integrating of binary orbits greatly reduced (Makino and 
Hut, 1990). The orbits of binaries with relative perturbation larger than 10~^ must be 
integrated numerically. To reduce the number of steps in integrating the orbits of such 
binaries, another technique — Kustaanheimo-Stiefel regularization — has been used. 

Over the last 20 years, the program NB0DY5 developed by Aarseth (1963, 1985, 
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1994) has been the most powerful tool for performing A'"-body simulations of star clusters. 
One important reason for the success of this program is that it can handle very efficiently 
the formation and evolution of binaries in A'^-body systems. A key part of NB0DY5 
is Kustaanheimo-Stiefel (KS) regularization (Kustaanheimo and Stiefel 1965), which is 
an extension of Levi-Civita's regularization of the planar Keplerian problem to three 
dimensions (Levi-Civita, 1956). In KS regularization, a Kepler orbit is transformed into 
a harmonic oscillator and the number of steps needed for the integration of an orbit is 
reduced significantly. 

In order to apply KS regularization, one first identifies a close pair of particles using 
some heuristic criterion, then integrates their center-of-mass and relative motion separately. 
The center-of-mass motion is integrated in the same way as the motion of all other particles 
in the system; the relative motion of the binary components is integrated using regularized 
time and coordinates. 

To reduce the accumulation of integration error, a technique known as energy 
stabihzation has been used (Baumgarte, 1973; Aarseth, 1985). The basic idea of energy 
stabilization is to introduce an artificial stabilizing term into the equation of motion so 
that the energy converges to the "true" value. In the case of an isolated binary, the energy 
is constant, and the force is simply adjusted so that the calculated energy converges to 
the original value. For a perturbed binary, the energy is integrated separately, and the 
force is adjusted so that the calculated energy converges to the integrated quantity. The 
stabilization technique vastly improves energy conservation; however, it does not conserve 
angular momentum. Thus, stabilization may give the wrong answer if the evolution of the 
binary eccentricity is important. 

Recently, numerical integrators which have no secular error in the integrals of motions 
have attracted the interests of both astrophysicists and numerical analysts. The integrators 
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which had been most extensively studied are the "symplectic" integrators (for an overview, 
see Sanz-Serna and Calvo, 1994), which have the property of being symplectic (or canonical) 
transformations in phase space. Since a canonical transformation is one way to define a 
Hamiltonian system, it is reasonable to suppose that an integration scheme which is itself a 
canonical transformation might describe the nature of the system better than an integrator 
which is not. In fact, symplectic integrators have the remarkable property that the errors in 
the integrals of the motion have no secular term. In the case of the Kepler problem, neither 
the energy nor the angular momentum change systematically in time. 

Time-symmetric integrators have not received as much attention as symplectic 
integrators. However, it has been demonstrated that some time-symmetric integrators also 
exhibit no secular errors, regardless of whether or not they are symplectic (Quinlan and 
Tremaine, 1990). One serious problem with both symplectic and time-symmetric schemes is 
that they work well only with constant stepsize. If variable stepsize is used, the performance 
of a symplectic integrator degrades to that of more "usual" schemes, in that it no longer 
automatically conserves the integrals of motion (Skeel and Gear, 1992; Sanz-Serna and 
Calvo, 1994). 

Recently Hut et al. (Hut, Makino and McMillan, 1995, hereafter HMM) have 
developed a novel technique to avoid the accumulation of the error in variable stepsize 
time-symmetric schemes. They found that the increase of the numerical error in the energy 
is dramatically reduced if the size of the timestep is determined in a time-symmetric way. 
When constant stepsize is used, a time-symmetric scheme is time-reversible, in the sense 
that, if we integrate the system forward for one timestep and then integrate it backward, 
the system returns precisely to the initial state. In the case of variable stepsize, however, 
the system generally does not return to the initial state, since the stepsizes used for the 
forward and backward steps are different. HMM used a time-symmetric criterion for the 



timestep, ensuring that the stepsize for the forward and backward steps are exactly the 
same. They apphed this scheme to the numerical integration of the Kepler problem in 
physical coordinates, and found that there was no secular error in the binding energy 

In this paper, we describe the implementation of this time-symmetric, variable step 
scheme to an A^-body code using the KS regularization. We have developed a fourth-order 
integration scheme and performed a series of numerical simulations. We find that our 
algorithm has the good properties of both symmetrized timesteps and KS regularization. 
First, like all symmetrized schemes, it shows no secular error in either energy or angular 
momentum. Second, like the code with KS regularization, it can handle binaries of arbitrary 
eccentricity with a constant number of timesteps per orbit. 

In section 2, we show the equations of motion and KS regularization. In section 3, we 
describe the integration algorithms used. In section 4, we describe the implementation in 
more detail. In section 5 we compare the results of our numerical experiments with those 
of other schemes, considering both isolated and perturbed binaries. As an example of a 
perturbed system, we consider a hierarchical three-body system in which the inner binary 
is integrated in KS coordinates and the outer binary is integrated in Cartesian physical 
coordinates. A summary and discussion are presented in section 6. 

2. Equations of Motion 

2.1. Cartesian Coordinates 

The equations of motions of the two particles to which we apply the KS transformation 
are expressed as 



where t is physical time and the gravitational constant G is set to unity. The 3-dimensional 
vectors rj and Vj are the positions and velocities of the particle i in Cartesian coordinates, 
rrii is the mass of particle i, and Fj is the force on particle i due to the rest of the system. 

In order to apply the KS transformation, we introduce the motion of the center-of-mass 
of two particles: 

miFi + m2r2 



mi + m2 
miVi + m2V2 

mi + m2 



(3) 
(4) 



and the relative motion: 



r 

dr 

'dt 



The accelerations are given by: 
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dt 



r2 - ri, 

V = V2 - Vi. 



(mi + m2)r 
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Fi + F2 

mi + m2 ' 



+ P, 



where the second term P in equation ([^ is the perturbation: 



Y = — - — 

1712 rrii 



(5) 
(6) 



(7) 



(9) 



2.2. KS regularization 



In KS regularization, the 3-dimensional position in Cartesian coordinates is transformed 
into a 4-dimensional position in KS coordinates. The physical time t is also transformed 



-8- 



into the KS time r. The transformation is expressed as follows: 



|r| ^ ul + ul + ul + u\, 



£(u)u, 



(10) 
(11) 



U2 Ui —U4 —Us 

y Us U4 ui U2 J 



(12) 



dt 

^ = r=(u,u), 



(13) 



where (u, u) is defined as the inner product of u and u. The resulting regularized equations 
of motion take the form 



u 



1 

2' 



" - ^/iu + ;r£^(u)P, 



1 

2' 



(14) 



where ' denotes differentiation with respect to r, i.e. ' = d/dr, superscript t denotes 
matrix transposition, and h is the specific binding energy of the two bodies, which may be 
expressed in the KS system as follows: 

2u' • u' — (mi + 1112) 



The time variation of h is described as: 



(15) 



h' = 2u' • C (u)P. 



(16) 



For an isolated binary, P = 0, so /i is constant. 



3. Integration Techniques 
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3.1. The Hermite Scheme 

The Hermite scheme (Makino, 1990) is a fourth order predictor-corrector scheme, 
which may be expressed as foUows: 

re = rfe + ^(ve + Vb)At-^(ae-ab)At2 + -i-(j^+j,)At3, (17) 

Ve = v, + ^(a, + a5)At- ^(je-j5)At2, (18) 

where ax and jx are the acceleration and the "jerk" (time-derivative of the acceleration), 
and subscripts h and e denote the values at the beginning and the end of the step, 
respectively. 

As can be seen from equations ([T7|) and (18), the Hermite scheme is an implicit scheme. 
The integration of one step is carried out as follows. 

[1] Predict the positions and velocities: 

Ypred = rfe + VfcAt + ia^At^ + ijj,At3, (19) 

vp.ed = Vfe + a;,At+ ijfcAt^. (20) 

[2] Evaluate the acceleration and the jerk je at the end of the step. 

[3] Correct positions and velocities using equations (|17D and (plSj) , using the following 
formula: 

re = r,.e. + ^r[^) At^ + ^rf A^^ (21) 
Ve = v,,e, + iri^)At3 + lrf A^^ (22) 

where r|,*^ are the i-th derivatives which are calculated as: 

rf^ = ^[-6(a6-a,)-At(4j, + 2j,)], (23) 
rf^ = ^[12(a6-ae) + 6At(jfe+je)] (24) 
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[4] Repeat [2]- [3] until the values of and Vg converge. 

Equations (pTD-(plD are equivalent to equations (|T^ and (p^. However, if the stepsize 
is calculated from higher order derivatives, equations (|21| ) -(p^ ) are preferred. 

The local truncation error of equation (|T^ is O(At^) and that of equation (^) is 
0{At^). The global error is 0{At^) for both r and v. One could also use the truncated 
form, 

re = rb + ^(ve + v,)At-^(ae-afe)(At)2, (25) 

which has a local error of O(At^) and a global error (that is, integrated over some fixed 
time interval) of O(At^) for both r and v. For the truncated Hermite scheme, the local 
truncated error in the position during one step is O(At^), while it is O(At^) in the full 
Hermite scheme. In either case, the global energy error is O(At^). 



3.2. Time-symmetrization 

We now describe the basic idea of the time-symmetric scheme. First, we consider an 
integration scheme with a constant stepsize. An arbitrary integration scheme with constant 
stepsize At may be expressed compactly as follows: 

ee = /(e6,ee,At), (26) 

where ^ represents the phase-space variables, i.e. C, = (r,v). If / does not depends on ^e, 
the scheme is explicit. Otherwise, the scheme is implicit. 

If the scheme is time- symmetric, the following equation holds 

e^/(ee,e;.-At)=6. (27) 
For example, a leap-frog scheme with constant stepsize, which is a well-known symplectic 
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scheme, is also time-symmetric. The leap-frog scheme is an explicit scheme. However, many 
useful time-symmetric schemes are implicit. 

The simplest example of an implicit time-symmetric scheme is the implicit trapezoidal 
rule, which may be expressed as: 

L = ^b+lmb)+m)w, (28) 

where 

0(0 - |. (29) 

In order to integrate for one step, the implicit equation (^) is solved iteratively until 
convergence is obtained. There are also time-symmetric integrators of higher order. The 
Hermite scheme [equations (|I^) and (|T^), or ( P^D and (|T^)] is an example of a fourth-order 
symmetric integrator. 

Now we consider time-symmetric schemes with variable stepsize. If the scheme is 
time-symmetric, the following equations are satisfied: 

= /(6,ee,Atb), ^^^^ 

6 = /fe,a,-Ate), 

where Atf, and Atg are the stepsize determined at the beginning of the step and the end of 
the step by some function S{C,), i.e. 

Atb = 5(6,ee); At, = 6{^e,^b). (31) 

In order to guarantee that equations ( PPD hold, it is sufficient to make the stepsize criterion 
time-symmetric : 

<^(6,ee) = 5(ee,6). (32) 



- 12 - 



One simple way to construct a timestep criterion that satisfies equation ( p2D is to take the 
average of stepsize calculated at the beginning and the end of the step. For example, 

S{C,,Q^^[si^,) + s{Q], (33) 

or 

S{^i.,Q^^l[s{^by + s{Q% (34) 

Equations ( p3| ) and (0) are implicit equations. The stepsize is determined by solving 
equation (|33| ) or (|3^ ) coupled with equations (|1^ and (0) [or ( pS] ) and (|T8|)]. 



3.3. Energy Stabilization 

Energy stabilization has been used in traditional programs to conserve the binding 
energy of a binary. In a stabilized scheme, energy conservation is enforced by introducing 
an artificial acceleration into the equation of motion, as follows (Baumgarte 1973): 

(fr _ (mi + m2)r ,-p_fOi\ {hnteg - hdirect)^ 

dt^~ |r|3 +^ [At) (vv) • ^^^^ 
Following the prescription by Aarseth (1985), equation (|35|) is transformed to KS coordinates 



as: 

72 



dr2 2 2 VAr/ M 

The third term of the right-hand side of equation (|36|) is the artificial stabilizing acceleration. 
The terms hinteg and hdirect are the specific binding energy of the binary obtained by 
integrating equation (|16]) and evaluating equation (|15|), respectively. Figure |l| shows the 
time variation of hdirect for different values of a. When a is less than 0.3, the behavior of 
the energy is simple. When a is larger than 0.3, the behavior of the energy is complex. 
With small a, hdirect cannot follow hinteg fast enough. In fact, for a = 0.1, a secular error of 
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~ 5 X 10~^/orbit was observed. For large a, the spikes around the periastron is very large 
and the error is enough adjusted every period. Aarseth recommended to use a — 0.4. In 
this paper, we set a to 0.4. 

4. Implementation of the Symmetrized KS Hermite Scheme 

4.1. The KS Hermite scheme 

The equations of motion of the KS binary are (14) and (16). We integrate them using 
the time-symmetrized Hermite scheme. 

ue = u, + l(uW + uW)Ar - l(u(^) - uf ))Ar^ + ^(uf ) + )Ar^ (37) 
u« = u« + i(uf +uf))Ar-l(uf -urVr^ (38) 
he = /i6 + ^(tf^ + M'^)Ar^^(/.f)-/if)Arl (39) 

Here and below we denote the i-th derivative with respect to a variable r by using the 
operator For example, 

(41) 

The second and third derivatives are calculated as (Aarseth 1985): 

uf) = hi^n,+ ]^nCl{u^)^^, (42) 

= l/.,u(^) + l/.«u, + \j-[nCl{n,)]P, + lr,£^(u,)r,J,,,, (43) 
/iW = -2u(^) • £'^(u)P, (44) 
= -2[u(2) . £^(u)P + u(^) • £'^(u(i))P + u(^) • >C^(u)rJp] (45) 
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The quantities Rh and Cb are, respectively, the separation of the binary in Cartesian 
coordinates and the transformation matrix at time r = r^. Because of the numerical 
difficulties associated with the singularity in equation (|15D as i? ^ (Aarseth, 1972), we 
integrate the binding energy h of the binary instead of calculating it from its definition at 
each step. 



4.2. Timestep Function 

The timestep symmetrization technique described in section 3.2 can be applied to an 
arbitrary algorithm for calculating the stepsize. 

The standard Aarseth formula (Aarseth 1985), which is being used in most iV-body 
programs, is expressed as 



At = s(u) 



|u{4)| 


|u(2) 


+ 


U(3)|2 


|u(5)| 


|u(3) 


+ 


U(4)|2 



(46) 



Here 7] is an accuracy parameter. The stepsize criterion is known to have good properties. 

For general A^-body problem, this criterion works very well. For a perturbed two-body 
problem with KS regularization, however, a more simpler schemes might be sufficient. We 
tried two formula. 

If we adopt the symmetrized scheme, we may be able to carry out an accurate and fast 
time integration without the standard formula. We can use simpler functions such as 

= Vi^^, (47) 



or 



At = s(u) 



as the stepsize function for the KS timestep At. 





U(2) 


U 


+ 


UW|2 


|u(3)| 


u{i) 


+ 


U{2)|2 



(48) 
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Here "simpler" means that equations (0) and (58) don't require the higher-order 
terms such as u^^^ and u*^^-*. It makes the convergence faster. 

Note that each of the functions (|l6|) and (^8]) gives a constant stepsize in KS coordinates 
for the time integration of a simple binary. For a time-symmetrized scheme, it is better to 
use the function which gives a nearly constant stepsize. This is because the estimation of 
the initial value of iterations becomes easier. 



4.3. Implementation in General A^-body Systems 

4-3.1. For a Few Body Problem 

We now describe the implementation of the symmetrized KS Hermite scheme in a 
full A^-body code. We consider here only a shared-timestep scheme, in which the whole 
system is advanced with a single global timestep. This method is sufficient for 3- or 4-body 
systems. 

The procedure for a single step of the time integration is as follows. At the beginning of 
the step, the relative position and velocity of the binary components in KS coordinates are 
known. The positions and velocities of other particles are given in Cartesian coordinates. 
The binary perturbation P and its time derivative Jp are first calculated in Cartesian 
coordinates. Then the relative position of the binary components is predicted in KS 
coordinates with KS time step At, while the center of the mass of binary and all other 
particles are predicted in Cartesian coordinates with physical time step At, which is a 
function of At and u. 

The stepsize At in the regularized system and the corresponding physical stepsize At 
must be calculated in such a way as to ensure the time- reversibility of the whole A^-body 
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system. We use the following formula to advance the physical time: 

ti = to + At, (49) 

At = T(Ar)=t«Ar + t?)^ + t(^)^. (50) 
^ ' 2 2 24 2 1920 ^ ' 

Equation (|50|) , which is a 5th order expression, is calculated by using the first through the 
fifth derivatives of the positions (i.e. all derivatives available during the integration). The 
formulation is the same as used by Aarseth (1985), except that we use the derivatives at 
the midpoint of the integration interval, which eliminates all even derivatives. Here t^f^ are 

2 

the i-th derivatives at r = Tq + 0.5Ar in KS time. They are calculated as: 



t^P = r = u ■ u, (51) 



where 



2 



t^P = 2(u-u(2) + u«-u«), (52) 

2 

tf = 2[u-u(^) + 4(u«-u(3))+3(u(2).u(2))], (53) 



ui - uo + — uo + — uo + + + , (54) 

u« = u« + — ^u^n ^uj^ (55) 

1 ° 2 ° 8 ° 48 ° 384 ° ^ ' 

(2) (2) , (3) , Ar^ 

uV = uV + — u^^ + -^u^^ + ^u^^ (56) 

uf = uin^u^n^u^^), (57) 

= ui^uf^^, (58) 



It may be sufficient to use the 4th-order Hermite scheme (the truncated Hermite 
scheme) to integrate the physical time, since the phase-space variables are globally 
integrated to 4th-order accuracy. More detailed error analysis will be given elsewhere. 
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The details of a single integration step are as follows. 
[1] Predict the positions and velocities of all particles. 

[2] Evaluate all accelerations and jerks using predicted positions and velocities. 

[3] Calculate the fourth and fifth derivatives of the relative position of the binary components 
in KS coordinates at both the beginning and the end of the step using the following 
formulae: 

(1) 1 
u, - 



- At 



6(uf)-uf) + 2Ar,(2uf)+uf)], (59) 

12rf)-uf)+6Ar,(ur) + uf)], (60) 

u(^) = ujnAr.uf, (61) 

u(^) = . (62) 

We calculate the higher order derivatives for other particles and the specific energy of the 
binary in the same way. 

[4] Calculate the new stepsize and the physical stepsize at the beginning of the step. For 
example, 

^Tnew = ^ [siUb) + s(Ue)] , (63) 

Atnew = T{ATnew), (64) 

[5] Correct the integrated values using the new stepsize Ar^eu, : 
Ue = U6 + uS'^Ar„e^ + iui')Ar2^^ + 

u(l) - U^^)+U^2)a ,1^(3) a 2 , 1u(4)a 3 , ^ (5)a 4 /gg^ 

[6] Repeat procedures [2] -[5] until both the stepsize and the integrated variables converge. 
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About the convergence of integrated variables, there is one simple method in which 
the procedure shall be repeated until the variables of all particles converge. Even if we 
adopt this simplest method, the number of iterations of the procedure is not so large for a 
few-body systems. In fact, in an example shown in the following section, the number of 
iterations was typically 4 for the case of an integration of a hierarchical triplet. The number 
docs not much exceed that required to converge the phase space variables in the ordinary 
Hermite scheme, which is 2 ~ 3. 

Even for more complicated cases, wc may make the number of iterations small if we 
choose an appropriate criterion formula to determine the stepsize At, or if we choose an 
appropriate accelerator to converge the variables. 

4-3.2. For a Large N Body Systems 

The outline of the procedure of integration of one step is same as that for the cases of 
few body systems. However, there are two difficulties. One is the increase of the number 
of iterations until convergence. Another is how to implement both symmetrized scheme 
and individual time step scheme. For the system with larger value of N, the number of 
iteration may become larger. For larger values of N, individual timesteps become necessary 
to achieve reasonable efficiency (Aarseth 1985). 

Implementation of time-symmetrization for individual timesteps will be discussed 
elsewhere. Here we only give a comment about this point. In a time integration of a 
large A'"-body systems, it is not necessary to symmetrize all variables of all particles. To 
symmetrize the integration of binaries (and its close neighbors) would be sufficient. 
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5. Numerical Experiments 

Figures 0a-^ show the behavior of the errors in the energy, eccentricity, and angular 
momentum of a binary with initial eccentricity e = 0.999999. In figures |^a-0c, the values at 
the apocenter are plotted. These figures show that there is no secular error for 2000 binary 
periods. The specific energy is calculated by equation (15). 

In the following subsections, we compare (1) the symmetrized KS-Hermite scheme, 
(2) the "plain" (i.e. neither symmetrized nor stabilized) KS-Hermite scheme, and (3) 
the stabilized KS-Hermite scheme, for e = 0.9. We have carried out other comparisons, 
with initial eccentricities ranging from e = 0.0 to e = 0.999999. In all cases we obtain 
qualitatively the same results as are shown below. 

We have investigated the dynamical evolution of both a simple binary and a hierarchical 
triple. The initial conditions used in our experiments are summarized in Table 1. In 
the case of the binary the integration was performed for 2000 orbits. For the triple, we 
integrated the system for 2000 periods of the inner binary. 

We used equation (ESf) as the stepsize function and equation (133) as the symmetric 



criterion. The coefficient rj in equation (|i8| ) was set to 0.01, corresponding to a stepsize 
of roughly 1/30 of the (inner) binary period. The condition to stop the iterations was 
|Are — Ar^l < 1.0 x 10~^^. The total number of timesteps was therefore ~ 6 x 10^ in all 
cases. All calculations were done in double precision (16-digit accuracy). 



5.1. A Simple Binary 



For this case, the number of iterations required by the symmetrization of stepsize Ar is 
2. The number of iterations necessary is small since Ar should be constant, i.e., Ar^ = Axe 
for the criterion (^). 
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Figures ^ and ^ show the result of an integration of the Kepler two-body problem 
for 2000 orbital periods of the binary. Figure ^ shows the error in the total energy (i.e. 
the error in the binding energy of the binary) for the three schemes listed above. Here 
we plot the error in the energy calculated from equation ([l5|) , not the error in the energy 
integrated using equation ([T6|) . The reason is that the integrated hinteg is constant, because 
the binary is unperturbed. The calculated energy thus indicates the accuracy of the orbital 
integration. We also use this calculated energy for other cases discussed below. 

Figure |^ shows that both the symmetrized and the stabilized schemes conserve energy, 
while the error in the energy increases linearly with time for the plain integrator. In the 
symmetrized case, the energy error AE/E is less than 10^^^ after 10'^ orbits, or about 10^ 
steps, indicating that it is mostly due to round-off error. 

We have plotted the change of energy at regular time intervals. The error of the 
energy becomes large at the pericenter, though the error returns back to nearly zero at the 
apocenter. As a result, the beating of the frequency of the binary orbit and that of our 
sampling is observed for the symmetrized and stabilized cases. The behavior of energy in 
the stabilized case in figure |^ is more complicated than that in the symmetrized case, since 
the time variation of error around the pericenter is complicated (see figure |I|d). 

Figure ^ is the same as Figure ^, but for the error in the angular momentum AA. 
Figure ^ shows that both the symmetrized and the plain schemes conserve total angular 
momentum, while the error in the angular momentum increases linearly in the stabilized 
case. The relative error in the angular momentum in the symmetrized case is again less 
than 10~i2_ 

Figure ^ is the same as Figure ^, but for the eccentricity. This figure makes explicit 
the fact that the eccentricity of the binary is conserved by the symmetrized integrator, but 
not by either the plain or the stabilized schemes. This result is quite natural since the 
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eccentricity depends on both the energy and the angular momentum. The plain scheme 
conserves angular momentum but not energy, while the stabilized scheme conserves energy 
but not angular momentum. Thus neither scheme preserves the eccentricity. Only the 
symmetrized scheme conserves both energy and eccentricity up to round-off error. 

5.2. A Hierarchical Triple 

For this case, the number of iterations required by the symmetrization is about 4. That 
required by the convergence of the phase space variables is 2 ~ 3 for the non-symmetrized 
Hermite scheme. The symmetrized scheme, therefore, is not much expensive than the 
non- symmetrized scheme. 

Figures ^ and |^ show, for our three integration schemes, the relative errors in the total 
energy IS.E jE (Figure ^ and the total angular momentum (Figure ^of the three-body 
system. In this case again, the relative error in the energy using the symmetrized scheme is 
no more than 10^^^ and the error in the angular momentum is about 10^^^ after 10^ steps. 

Figures |8|-|TT| show the evolution of the inner binary. From the analytical treatment 
presented in the Appendix, the time evolution of energy, angular momentum and eccentricity 
of the (weakly perturbed) inner binary may be expressed as follows. 

^EjE = 0.0, (67) 
AA = 0.5 ■ 10"^[1 -cos(2not)], (68) 
Ae ^ 0.25- 10"^[-1 + cos(2fiot)], (69) 

where subscript indicates unperturbed values. 

Figure ^ shows the variation of the specific energy of the binary. The specific energy 
is conserved in both the symmetrized and the stabilized cases, while it increases linearly in 
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the plain case. Figure ^ is the same as Figure ^, but for the binary angular momentum. 
Figure |T0| is the same as Figure |, but for the binary eccentricity. In figures |^ and [1^, 



the bare values at the apocenter of the inner binary are plotted. As shown in figures |8HlO|, 
since the amplitude of the variations during one period of the outer binary is much larger 
than the numerical errors incurred within one period, as shown in equations (64)- (66), the 
numerical error could not be seen if the bare values were plotted. In order to show the 
secular variation, we applied the linear regression. Straight lines in figures P|-pil| are the 
results of the least square fitting. Figures ^, ^ and 10 show that there are no significant 



secular variations in the symmetrized case, and that the eccentricity of the binary is not 
correctly followed by the stabilizing integrator. Comparing figures |^ and ^ we see that the 
error in the total angular momentum comes mostly from the integration of the KS binary 
in the stabilized case. 



Figures [Tl| a and b show the time variation of the angular momentum and the 
eccentricity of the inner binary during the first few periods of the outer binary. The dotted 
curves correspond to the analytical values, while the solid curve corresponds to the results of 



numerical experiments in the symmetrized case. Figures ^Tp. and b show that the numerical 
experiments agree well with the analytical results. 

To summarize, our numerical experiments show that both the total system and 
the inner binary are integrated correctly in the symmetrized case. In other words, the 
momentum transfer between the outer particle and the inner binary is accurately followed. 
Thus, our scheme follows the total system consistently, even though part of the system is 
expressed in KS coordinates and the rest in physical coordinates. 
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6. Summciry and Discussion 

We have developed a new algorithm for the integration of perturbed binary motion. The 
algorithm is constructed using two techniques: KS regularization and time-symmetrization. 
Our algorithm allows variable stepsize, but produces no secular error in either the binding 
energy or the eccentricity We have presented the results of numerical experiments for one 
case with e = 0.9. We have also performed experiments for other eccentricities, and have 
confirmed the superiority of the method in cases ranging from circular (e = 0) to highly 
elongated (e = 0.999999). 

Our scheme is composite, in the sense that the internal motion of the binary and the 
motion of the rest of the system are integrated in separate coordinate systems. This is 
rather different from the usual applications of symmetric or symplectic integrators, where 
all degrees of freedom are integrated in the same way. For an unperturbed binary, our 
algorithm simply applies the symmetrization to the system with a transformed Hamiltonian. 
Thus it is not surprising that our algorithm works well for unperturbed binaries. 

For perturbed systems, the result of numerical experiment shows that the error caused 
by any inconsistency between the integration of the system in the physical coordinates 
and the KS coordinates is smaller than the error caused by the integration of KS binary. 
Furthermore, our numerical experiments shows that the calculation cost of the new scheme 
is not very high compared with that of ordinary Hermite scheme. This result suggests that 
our scheme is valid for weakly perturbed cases. An accurate integration of long surviving 
binary is the main target of the time symmetrized scheme. Since the perturbation on such 
a binary is weak, our scheme is excellent from a practical point of view. 

In conclusion, we have shown that the time-symmetrization scheme, introduced by Hut 
et al. (1995) can be successfully generalized to include KS regularization. It is clear that 
our scheme can provide excellent accuracy for individual demanding cases. Under which 
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circumstances this new scheme will prove to be competitive will have to await a detailed 
testing in realistic simulations. 
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partly by Joint Research Program Between the Royal Society, the British Council and 
JSPS. 



A. Appendix: Analytical treatment of a hierarchical triple 

We introduce physical coordinates as shown in figure 12. The origin is placed at the 
center of mass of the inner binary. The masses of particles 1 and 2 are mi = m2 — m — 1/2. 
The unperturbed motion of particles 1 and 2 may be expressed as 

ri = (ri cos 0, rising), (Al) 

Vi = (fi cos — ri0sin 0, fi sin + ri0cos0), (A2) 

r-2 = -(r2COS0,r2sin0), (A3) 

V2 = — (r2 COS0 — r20sin r2 sin + r20cos0), (A4) 

where r, and Vj are the position and velocity and the radius of particle i, — |rj|, and (f) is 
the true anomaly of the inner binary. 

The unperturbed motion of particle 3 is expressed as a circular orbit around the center 
of mass of the binary: 

ra = [R cos ip^Rsinip) = [R cos Vtt, R sin Vtt) , (A5) 
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where R is the distance from the third body to the binary center of mass and ijj is the true 
anomaly of the outer orbit. 

Since mi = m2, ri — r2 — r, where r is half the distance between particles 1 and 2. 

The accelerations of particles 1 and 2 due to particle 3 are: 

(Rcosib — r cosd) Rsinib — rsind)) 

ai = 777.3^ 3 — -, (A6) 

[i?2 + r2-2i?rcos((/.-V')]i 

(i?cos'0 + r cos f/), i?sin-0 + r sin 0) , . „, 

a2 = rris- — ^— ^. (A7) 

[R^ + + 2Rr cos{(j) - ij)]^ 



2 



For unperturbed motion of the inner binary, the binary separation and time derivatives 
of the separation and the phase angle are 

2(l + ecos(/))' ^ ^ 

Ae sin 

- ^ r ^ ^^(13^' (^9) 
. A A{l + ecos<Pf 

Here a, e and A are the semi-major axis, the eccentricity, and the specific angular 
momentum of the relative motion of the inner binary; A = \Ja{l — e^). 



A.l. Energy Variation 

The specific binding energy of the inner binary, and its variation due to infinitesimal 

changes 5xj and Svi, are 

E - E^-v,-^, (All) 

i=l,2 ^ ^' 

5E = 5^v,-5v, + ^^^ = -$:v.-'^v„ (A12) 

i=l,2 i=l,2 

Expressing the variation of the velocity as 5vj = SiiSt, we obtain 

6E = - E ■ ikSt. (A13) 

j=l,2 
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From equations (|A1|) -(|7^)), we obtain the following equation: 



dE_ 



Rr cos(0 — ■?/') — Rrcj) sin(0 — ip) — rr 



[R^ + - 2Rr cos{(f) - 
Rr cos{(f) — tp) — Rrcj) sin(0 — ip) + rr 



[R^ + r2 + 2Rr cos( 



(A14) 



Substituting equations 
respect to a/R, we obtain 



and ( |A10| ) into equation ( |A14|) , and expanding with 



dE f . 6r cos(0 — ip) 
= < — 27^?^ -I- 

dt R^\ R 



Rr cos {(p — ip) — Rrcp sm{(p — ip) > . (A15) 



Assuming that the change in the position of particle 3 is negligible, we can integrate 
equation ( [A15| ) along the unperturbed orbit of the inner binary. Integrating it over one 
period of the inner binary, we obtain: 



AE 



0. 



(A16) 



This result is usually called as adiabatic invariance. 



A. 2. Angular Momentum and Eccentricity 



The change in the specific angular momentum of the inner binary is 



A 
SA 



2 ^ |ri X Vil, 



=1,2 



1=1,2 



(A17) 
(A18) 



Again we neglect the variation in R as we did in the previous section and replace 6v with 
a6t, to find 



dA = 2 



|a(l — e^) [ [(cos (p, sm(p) x (Rcos ip — r cos (p^Rsmip — r sin > 



1 + e cos ( 



[i?2 + r2 - 2Rr cos( 
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[— (cos sin 0) X [R cos i/) + r cos (j), R sin + r sin ( 



2 mg- 



[i?2 + r2 + 2i?r cos( 
ii?a(l - e2) 



sin('?/' — 0).?? 



;i + e cos (f)) [R^ + r2 - 2i?r cos(0 - i 
Expanding equation ?? with respect to a/i?, and neglecting higher terms, we obtain: 



dA 
lit 



2 ms- 



a(l-e' 



2/?2 



sin('?/' — 0) 



1 + e cos ( 



6a(l — e2) cos(0 — ■?/') 



2i? 



1 + e cos ( 



(A19) 



(A20) 



Assuming that the position of particle 3 does not change during one period of the 
inner binary (i.e. R and ip are constant), we integrate equation ?? over one period of the 
unperturbed motion of the inner binary: 

3ms a^(l — e2)^ 



AA 



i?3 



simjj cos ipB(e). 



(A21) 



From the variation of energy and angular momentum, we obtain the variation of 
eccentricity for one period as follows: 



Ae 



:i-e2)5 ^ ^ a(l-e2) ^ ^ 
^ r^AA + ^-AE 



where the function Bie) is defined as 



Bie) 



ea2 

3 



{3(1 -e^)^5(e)}sin7/'cosV', 



2^^ 2 cos2 6-1 



lo (1 + ecos0)^ 
For our numerical experiments, the value of B{e) is -6(0.9) = 4260. 



(A22) 



(A23) 



A. 3. Secular Variation 



Integrating equations ( |A16 ), (|A21j ) and ( A22 ), the variations of the specific angular 
momentum and the eccentricity of the inner binary are 

= 0.0. (A24) 



E 
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7 , 



AA{t) = ^ f ^) ^Iii-^5(e) [1 - cos(2fiot)] , (A25) 



Pin / 

3n ^2\4 



Here Pj^ is the unperturbed period of the inner binary. 

For the case of our numerical experiments (i.e. e = 0.9), equations ([A24| )-( |A26|) yield 

^ - 0.0, (A27) 
AA{t) ^ 0.46 • 10^^(1 - cos 2^0^), (A28) 
Ae(t) ^ 0.23- 10"^(-1 + cos 2^]o^). (A29) 
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Table 1: Initial orbit parameters for the simple binary case. 



mi 1712 ah e period[27r] 

0.5 0.5 1.0 -0.5 0.9e+0 1 



Table 2: Initial orbit parameters for the hierarchical triplet case. 





"'1 


1)1-2 


(I 


// 


(■ 


])(Ti()(l[27r] 


inner 


0.5 


0.5 


1.0 


-0.5 


0.9e+0 


1 


outer 


1.0 


0.01 


10.1 


-0.05 


O.Oe+0 


32 
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Fig. 1. — Effect of stabilizing parameter a on time variation of error in binding energy of a 
binary with eccentricity e = 0.9. Figures (a) through (f) correspond to the cases of o; = 0.1, 
0.2, 0.3, 0.4, 0.5, and 0.6, respectively. 

Fig. 2. — Time evolution of errors for the case e = 0.999999. Unit of time is the period of 
the binary. 2a). Relative error in the specific energy of the binary. 2b). Error in the angular 
momentum. 2c). Error in the eccentricity. 

Fig. 3. — Time evolution of relative error in binding error for the simple binary case. Unit 
of the horizontal axis is the period of the binary in KS coordinates. Solid, dashed and long- 
dashed lines correspond to the cases of symmetrized, "plain" and stabilized case, respectively. 

Fig. 4. — Same as Figure 3 but for the angular momentum. Top and bottom are the same 
figures except for the scale of vertical axis. 

Fig. 5. — Same as Figure 3 but for the eccentricity. 

Fig. 6. — Time evolution of the error of total energy for the perturbed binary case. The 
curves are same as those in Figure 3 The unit of horizontal axis are the period of the inner 
binary in its KS coordinates. 

Fig. 7. — Same as Figure 6 but for the angular momentum. 

Fig. 8. — Time evolution in the specific energy of the inner binary for the perturbed binary 
case. The curves and and the unit of horizontal axis are same as those in Figure 6. 

Fig. 9. — Same as Figure 8 but for the specific angular momentum of the inner binary. 

Fig. 10. — Same as Figure 8 but for the eccentricity of the inner binary. 

Fig. 11. — The evolution of the inner binary for the symmetrized case. Unit of time is the 
period of the inner binary. Solid and dashed curves correspond to the numerical experiment 
and theoretical, respectively. 11a). Angular momentum, lib). Eccentricity. 



